Quantum optomechanics beyond linearization 
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The quantum dynamics of optomechanical systems was mostly studied for their fluctuations 
around classical steady states. We present a theoretical approach to determining the system ob- 
servables of optomechanical systems as genuine quantum objects, for example, a coupled quantum 
mechanical oscillator to a cavity single photon. In this approach we study the dynamics of such 
systems in strong coupling regime. We find that, under strong optomechanical coupling, steady 
quantum states of optomechanical systems driven by continuous-wave single photons exhibit peri- 
odic oscillation and cavity noise considerably affects system observables. 
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I. INTRODUCTION 



Optomechanical systems (OMS) provide a platform to 
study macroscopic quantum phenomena [l|, l2| . Recent 
experiments have achieved the motional ground states of 
micro-nano oscillators [3-5], creating the possibility to- 
ward manipulating macroscopic quantum objects. On 
the other hand, experimental systems are approaching 
the single-photon strong coupling regime [6j-l8|, where 
the radiation pressure of a single photon could displace 
mechanical oscillator by more than its zero-point un- 
certainty. Strong optomechanical coupling is useful to 
making the macroscopic superposition states proposed 
for testing quantum theory |9l-[llj. 

The theoretical understanding on this novel regime 
of single photon and strong coupling optomechanics 
is under way too. The initial studies by Rabl |13 |. 
Nunnenkamp and co-workers [l3j, apply the linearized 
Langevin equation and/or the master equation about 
quantum fluctuations around classical steady states to 
find system observables. Another existing theoretical ap- 
proach adopts the Schrodinger equation under the pure 
quantum state assumption which is valid in single-photon 
subspace and for negligible loss of mechanical oscillator 
[lj, [Da]. Other study involving strong optomechanical 
coupling can be found in [161 ]. Despite the progress, more 
appropriate treatment should be developed for quantum 
OMS that can not be simply described as classical steady 
state plus quantum fluctuation to have linearized equa- 
tions of motion. 

To have a clearer picture of the problem, we refer to the 
approach in the previous studies. This is the linearization 
of the Langevin equation from the Hamiltonian (K = 1) 
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of a generic OMS, where g is the optomechanical coupling 
intensity, E the driver intensity, and u) c , ui m the cavity 
resonance frequency, mechanical frequency, respectively. 
Here we consider a continuous- wave driver with the cen- 
tral frequency u>q. In a frame rotating with the driver 
frequency ujq and diagonalizing the system Hamiltonian 
(the terms except for the driving term in (U}), the exact 



Langevin equation for the cavity mode reads 



•— a + lAa + i a'aa 

2 UJrrt. 



■e-^iE-Mc), (2) 



where A = g 2 /uj m — Ao with Ao = uj c — ujq being the 
detuning of the driver frequency from the cavity fre- 
quency, Q m = ig/ui m (b' — b), and k, £ c the cavity damp- 
ing rate and cavity noise operator, respectively. The ef- 
fect of the nonlinear term, the third term on the right- 
hand side of dSl, is canceled with the argument of weak 
driver in [l2l Il3| . In fact it is to perform a displacement 
a — » a s + a, where the classical steady state value a s 
of the cavity field is proportional to the driver intensity 
E, so that Eq. @ will be reduced to a tractable linear 
equation neglecting the higher order terms of E. The 
results based on such linearization are consistent with 
the numerical calculation with the master equation about 
the fluctuations around classical cavity steady states [l3| , 
showing its equivalence to the practice of decomposing 
the system operators into classical steady state values 
plus quantum fluctuations (see, e.g. [l8j ) adopted in 
most previous studies on OMS. However, in the situa- 
tions of a cavity single photon coupled to a mechanical 
oscillator in quantum regime, there is no way to define 
a classical steady state for the single quanta a since it is 
a quantum fluctuation itself. Then, for the single pho- 
ton optomechanics, the effect of the nonlinear term in 
@ should not be simply neglected to the dynamics of 
OMS, and any linearization for this type of problems is 
not appropriate. 

In this work we present an approach that models OMS 
as quantized interacting oscillators (one cavity and one 
mechanical mode) which damp energy into their envi- 
ronment and are subject to quantum noise effects at the 
same time. This full quantum approach allows one to 
determine the system dynamics without the reliance on 
the classical configurations of OMS, because it is based 
on the evolution of quantum operators. Though the ap- 
proach applies to the general coupling of OMS, we will 
primarily study the physics of quantum OMS in strong 
coupling regime by illustrating the evolution of the rele- 
vant system observables and clarifying the physical fac- 
tors to determine their values. Moreover, the effects of 
quantum noise can be well captured in the approach. 



II. SYSTEM EVOLUTION AND OBSERVABLES 

We start with the quantum state evolution of a generic 
OMS governed by the master equation (ft = 1) 

p = -i[Hs(t),p] + KV[a]p + -y m {n th + l)V[b]p 

+ i m n t hV[P]p = L{t)p (3) 

with the system Hamiltonian [171 ] 

H s (t) = -^ga)a{x m cos(uj m t) + p m sm(u rn t)} 

+ iE(e lAot a<< - e- lAot a) (4) 

in the interaction picture (the purpose for adopting the 
interaction picture will be seen below), where n^ — 
^ e u m /k B T _ -|^-i j g ^g thermal phonon number at the 
temperature T, and k, j m the cavity and mechani- 
cal damping rate, respectively. Here we define x m — 
(6 + 6^)/v / 2 andp m = — i(b— w)/V2 as the dimensionless 
displacement and momentum operator of the mechanical 
oscillator with the frequency cu m . T>[c)p = cp& — (pc ! c + 
c'cp)/2 is the dissipator in Lindblad form. The initial 
state of OMS is prepared as p(0) = \0) c (0\® Pth, the ten- 
sor product of cavity vacuum and mechanical thermal 
state. As in solving the Langevin equation in ([2]), the 
nonlinear coupling term in (U) makes it difficult to find 
the analytical solution to the master equation in ([3]). 

Here we provide a method to calculate the expecta- 
tion value (Os) = Trs{Osp(t)} of a system operator Os 
without the solution p(t) to the master equation in ([3]). 
To do this, we look at the combinations of OMS and their 
reservoirs. The damping of cavity field and mechanical 
oscillator can be described in terms of a linear coupling 
of the system operators Ci = a, 6 with the quantum noise 
operators £ c , £ m of the cavity and mechanical reservoir, 
respectively. In this picture the overall evolution for the 
combination of an OMS and its associated reservoirs is 
determined by the following unitary evolution operator 
of stochastic Hamiltonian [19| 



U(t,0) = Texp{-i / H s {T)di 
Jo 



+ Y.V7i (dBl^Ci-dBiirH)}, (5) 
Jo 

where T stands for a time-ordered operation, Bi(t) = 
Jo £i( T )^ r J li = K i7m! and Hs(t) takes the form in the 
interaction picture. The stochastic operators dBi and 
their conjugates satisfy the Ito's rules corresponding to 
the correlators (i c (t)H(t')) R = 6(t - t\<£(*)&(f)> R = 
of the cavity vacuum noise, and (tm(t)tm(t')) R = 
(n th + l)S(t-t'), (il(t)£ m (t')) R = nthS(t-t') of the me- 
chanical thermal bath; see Ref. [l9| for the details. Trac- 
ing out the reservoir degrees of freedom in the increment 
by infinitesimal transformation, U(t + dt : t)p(t)U^ (t + 
dt, t) — pit), restores the master equation in ([3]). The den- 
sity matrix p(t) for the combination of OMS and reservoir 



is assumed to factorize at t = 0, i.e. /5(0) = p(0) ® pr 
with pr being the tensor product of cavity reservoir vac- 
uum state and mechanical reservoir thermal state. The 
expectation value of a system operator can be rewritten 
with U(t,0) as 

Tr s {6 sP {t)} 
= Tr s {d s Uo(t,0)Tr R (U(t,0)p(0) ® p R U^(t,0))U %0)} 

= TrsMuHt,0)U o 1 (t 1 0)dsUo(t 1 0)U(t,0)p(0) ® PR }, (6) 

where Uo(t, 0) — exp{— i(uj c d^a + u> m wb)t} converts the 
state in the interaction picture to that of the Schrodinger 
picture. 

For any operator 0$ = /(ci,q), the transformation 
by U o (t,0) only adds the phase e~ iulit (e iUit ), where 
i = c or m, to Ci (q). Its expectation value as a sys- 
tem observable is therefore the average of the trans- 
formed functional operator W(t,0)f(ci,cl)U(t,0) = 
f(U^(t,0)ciU(t,0),U^{t, 0)c\U(t,0)) over the total ini- 
tial state of system plus reservoir, with the above- 
mentioned phase absorbed in Cj (cj). The reduction of de- 
termining system observables to finding the transformed 
basic system operators U'(t, 0)ciU(t, 0) is the main ad- 
vantage of our combined unitary evolution approach. 
The more complicated averages such as (Oi(£ + t)02(t)), 
which is gen erally calculated with the quantum regression 
formula [20J|, can be calculated in a similar way. 



III. DECOMPOSITION OF 
SYSTEM-RESERVOIR EVOLUTION 



We will find the expectation value of a system op- 
erator Os — f{ci,c\) by averaging its transformation 
/(C/ t (t,0)c l L>(t,0),L> t (t,0)c l t L>(i,0)) over the initial sys- 
tem state plus reservoir state. The total unitary evo- 
lution operator U(t,0), however, involves the driving 
on cavity, the optomechanical coupling, as well as the 
coupling between system and reservoirs, which are not 
mutually commutative quantum processes. To simplify 
the calculation of system operator transformations under 
such overall unitary evolution, we will need the following 
decompositions 

Te - l / t dr(ff 1 (r)+_ff 2 (r)) 

where V 2 (t,r) = Texp{-i/ r ' dr'-fl^r')}, and 

T e -i So <ir(Hi(T)+H 2 (r)) 

= Te-'itidTHiir) r Te -iftdTV?(T,0)H 2 (T)V 1 (Tfi)^ /g) 

where V\(t, 0) = Texp{— i L dT'Hi(r')}, of a time- 
ordered exponential. The proof of the two decomposi- 
tions is given in Appendix A. 



First, applying (J7J to ([5]), we separate out 
a unitary evolution operator Vd(£, 0) of the 
system-reservoir coupling so that the overall 
unitary evolution operator can be decomposed 
into the form U(t,0) = V s (t,0)V D {t,0), where 

V D (t,0) = Texp{E,(jfo y/TidBfrfa-fi y/TidMM)} 

describes the coupling between the system and 
reservoirs. The unitary operation VD(t,r) inside 
V S (t,0) = Texp{-if*dTV D (t,T)H s (T)Vl(t,T)} trans- 
forms the system operators c% — a, b in Hs(t) to 

V D (t,T)CiVl(t,T) = e-^-^Ci + hifr^^Cifar), 

(9) 

with hi(t,T) = y fTiJr e ~' ldt '~ T)/ ' 2 &( t ') dt ' bein S the in- 
duced quantum colored noise operators satisfying the 
commutation relation 



r 4 (T,r') = [hi(t,r),nl(t,T')] 



-7i|r-r'|/2 



e —n(t-r)/2 e —n(t-r')/2. 



(10) 



see Appendix B for the proof. From now on, one will 
rewrite the expectation value of a system operator as 

Tr s (6 sP (t)) = Tr StR {v£(t,0)Ut(t,0)d s U (t,0)Vs(t,0) 
x V D (t,0) P (0)^ PR V^(t,0)}. (11) 

This expectation value will be therefore determined 
by the transformations Vg(t, 0)ciVs(£,0), together with 



Vl(t,0)ciV D (t,0) that are similar to Eq. © [21 1 



Next, using Eq. ([8]), we will separate the pure driving 
process and the optomechanical process in the unitary 
evolution operator Vs(t, 0) = Texp{— i J Q d,TH' s (t 7 T)}, 
where H' s (t,r) = iE{a^(t,T)e lA ° T - a(t,T)e- lA ° T } ~ 
gK m (t,T)a^ (t,r)a(t,T), to have the product 
V s (t,0) = V E (t,0)V OM (t,0). The pure driving op- 
erator, V E (t,0) = Texp{J5/ t dr(ote iAoT e-§( t - T 5 - 
h.c.)} x Texp{£;/ *dT(e lAor nt(i,T) - h.c.)}, consists 
of the actions by the external driver alone and the 
noise h c {b, r) from the vacuum reservoir. The form 
of the optomechanical coupling evolution operator 
F OM (*,0) = Texp{igJ*dTK m (t,T)CHt,T)C(t,T)} 
is obtained by (JTJ) and ©, from which we have 
K m (t,T) = cos(uj m T)X m (t,T) +sin(w m T)P m (i,r) with 

X m (t,r) = V2e- 1 ^^- T) x m + {n m (t,T) + h.c.) 
P m (t,r) = V2e- 2 ^^ t -^p m - (in m (t,T) - h.c.) (12) 

being the transformed displacement and momentum op- 
erator of the mechanical oscillator, as well as the trans- 
formed cavity operator 

C(t,r) = vUT,0)V D (t,T)aVl(t,T)V E (r,0) 



= P -f(*-r) 



(a + D 1 (T))+n c {t,T)+D 2 {r) 



; {t ~^a + G(t,T). 



(13) 



In (Tl3l) . the functions due to the displacement by the 
external driver and cavity noise operator in V E (t, 0) are 
respectively found as 



D X { T ) = ( e -«(*-T)/2+<A r _ e -Kt/2) 

D 2 (t) = E f dt' e jAot 'r c (i',r). 
Jo 



E 
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The unitary evolution operator for the combination 
of OMS and reservoirs has now been decomposed as 
U(t,0) = V E (t,0)V O M(t,0)V D (t,0). Such decomposition 
of a quantum physical process into three dependent but 
simplified ones makes it possible to find the transforma- 
tions U^(t, 0)ciU(t,0) of the basic OMS operators. As 
it will be shown below, the contributions to a system 
observable from the different factors, e.g., the external 
driving and the optomechanical coupling, can be seen by 
such decomposition as well. It is therefore convenient for 
the approach to study the system dynamics in various 
different regimes, as one adjusts the system parameters 
to the decomposed processes. 



IV. DYNAMICAL AND NOISE EFFECTS 

Now we go back to study the physics of quantum 
OMS with the above decomposition technique. The in- 
finite product expression for the unitary evolution op- 
erator V O M(t,0) in V S (t,0) = V E (t,0)V OM (t,0) enables 
one to obtain an analytical form of the transformation 
Vg(t,0)dVs(t,0) (see Appendix C for the details). For a 
weak drive of E/k <C 1, its average over the initial cavity 
vacuum state reads 




400 



FIG. 1: (color online) Proportion of mechanical noise cor- 
rected cavity photon number R m = Anp/np v.s. quality 
factor Q = u) m / r Ym. The thick solid (orange) line is for the 
resolved-sideband parameter u) m /i^ — 1, and the thin solid 
(blue) line for io m /n — 2. The parameters are g/n — 0.5, 
Ao/w m = 0, and T = 0. The plots show the mechanical noise 
correction at nt = 10. 



(Vkt,0)aV s (t,0)) c = D x (t)+ig dr exp{-ig 2 8(T)} exp{i 5 / du e~ K(t ~ u) sin(o; m M) J P m (i, «)} 



x exp{z ff / du e- K{t ~ u) cos(uj m u)X m (t,u)}k m (t,T)e- Kit - T)/2 Gc{t,T), 
Jo 



2 



(15) 



where 

G c (t,T) = G(t,r) + ig I dTir c (Ti,T)K m (t,Ti)G(t,Ti) + (ig)* I d-nT c (Ti,T)K m (t,Ti) I dT 2 T c (T 2 ,T 1 )K m (t,T 2 )G(t,T 2 ) 

J t J r J n 

+ (ig) 3 /" drir c (ri,r)^ m (t,ri) / dT 2 r c (T 2 ,n)K m (t,T 2 ) f dT 3 T c {T 3 ,T 2 )K m (t,T 3 )G{t,T 3 ) + ■ ■ ■ (16) 

I 



includes the correction to the drive opera- 
tor G(t, t) defined in ([T^]) by the induced 
cavity colored noise h c with the correlation 
function T c (t, t'). The extra phase S(r) = 

2 J Q T due - ^* - ") sin(w m u) /„" dve"^ 1 ^^^^-^' 2 cos(uj m v) 
is due to the non-commutativity between X m (t,r) and 



X, 



P to (*,t). The term L>i(t) =£;(e lAo *-e- K */ 2 )/(f+iAn) 
simply arises from the pure driving process Ve(£)0). 

A. Weak coupling limit 

In the weak coupling limit g <C k, the average 
(Vg(t, 0)aVs(t, 0)) c could be approximated by the term 

„E * n pi(A +w m )t 

U ~ f + iAo 1 ^ ^ J MAo + ^) + (f + f) 



i(Ao-u m )t 



z(A 



+ (! + ¥)' 



(17) 



in the lowest g order of Eq. (|15p . in addition to the pure 
driving and noise correction terms. Showing the anti- 
Stokes and Stokes sidebands, this approximation well 
explains the physics of OMS in linearized and weak- 
coupling regime. 

An important feature in our approach is that a system 
observable should be found as the average of a trans- 
formed system operator over the initial state of both 
OMS and reservoir. The average over reservoir state 
manifests the mechanical and cavity noise corrections to 
the system observable. For instance, the colored mechan- 
ical noise operator h m in X m (t, r) and P m (t, r) of (fT5T> 
corrects the mean cavity photon number rip = (a^a). In 
the lowest order, which well describes the weak couping 
regime, it changes the photon number by 



Anp = 



9 2 E 2 
$ +A 2 Jo 



dn 



d T2 e-«(*-^)/ 2 e -"(*-^)/ 2 



x cos Wm (r 1 ~r 2 )e lAori - tAoT2 r m (T 1 ,r 2 ) (18) 

at T = (extra terms containing nth wm be added 
for a system at the temperature T > 0). The correc- 
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FIG. 2: (color online) Real-time evolution of cavity quadra- 
ture X c with the dimensionless time nt. The thick solid (or- 
ange) line shows the process for the coupling gjn = 0.1, the 
dashed line for g/n — 0.5, and the thin solid (blue) line 
for g/n = 1. Here we choose the system parameters to be 
E/n = 0.01, oj-m/n — 1, cu m / r y m — 100, A /cj m = 0, and 
T = 0. 



tions in the higher orders can be calculated with Wick's 
theorem to sum up the products of (n^(t, T)n m {t, t'))r 
and (n m (t, r)nj 7l (i, t'))r. The contribution of these cor- 
relators to photon number np is rather small under 
the condition "f m <C K. Fig. 1 illustrates the propor- 
tion of such correction, R m = Anp/np where np = 
((U+Dl)(l+Di))+Anp, in the total mean cavity photon 
number of the weak coupling regime, showing that Anp 
will become more negligible with the increasing quality 
factor of mechanical oscillator. 



B. Transition from weak coupling to strong 
coupling regime 



Before we discuss the cavity noise effect, we take a 
look at the effect of optomechanical coupling on cavity 
observables as an OMS undergoes the transition from 
weak coupling to strong coupling regime. Here we take 
the cavity quadrature X c = {a + a^)/v2 for illustration. 
From Eq. (|15p the further averages are taken over the 



initial mechanical oscillator and reservoir state to obtain 
the numerical values of X c . Fig. 2 shows the real-time 
evolution of X c for three different optomechanical cou- 
pling intensities. For the weakest coupling, the pure 
drive process Ve(£, 0) is dominant in the total process 
V s (t,0) = V E (t,0)Vo M (t,0), and X c will finally tend to 
\f2Re{as\, where a s corresponds to the classical steady 
state value for a. In the limit 5 = 0, the system will reach 
a static steady quantum state, with the cavity being in 
a coherent state, after a sufficiently long time. More sig- 
nificant effect of Vom(£,0) with increasing coupling in- 
troduces the periodic oscillation patterns to the function 
X c (t), after this observable becomes stable. Given dif- 
ferent detuning of the external driver, the stably oscillat- 
ing quantum states p s (t) of OMS differ greatly in strong 
coupling regime. This can be shown by the quadrature 
functions X c (t) = Tr s {p s (t)(a + a+)/V^} in Fig. 3. 




FIG. 3: (color online) Real-time evolution of cavity quadra- 
ture X c in strong coupling regime. The thin solid (blue) line 
shows the oscillation of cavity quadrature for Ao/um — 0; 
the thick solid (orange) line for the detuning Aa/uj m = 1; 
the dashed line for the detuning Ao/u; m = — 1. Here the 
system parameters are E/k = 0.01, g/n — 2, oj m /« = 1, 
tomhm = 100, and T = 0. 

The specific superposition of the harmonic components 
rujj m (the first harmonic appears in Eq. (|17p) makes the 
shape of a periodic pattern in Figs. 2 and 3, as its over- 
all repetition period is controlled by the mechanical fre- 
quency uj m . This feature is unique to quantum OMS in 
strong coupling regime, and reflects the fact that stable 
quantum states of OMS are not static, i.e. p(t) ^ 0, in 
the regime. Here we provide an interpretation for this re- 
sult with the OMS state evolution neglecting damping to 
environment {k and 7 m are very small to be neglected). 
For a zero temperature OMS driven at the detuning 
A = 0, its pure quantum state Te~ 4 Jo Hs ^ dT \0) Ct7n un- 
der such condition can be approximated by the state 

t 
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E(e 



i<t> m {t) 



f e-^ m(T) e- 4Xm(t ' r W)a}|0) c 
Jo 



(19) 



FIG. 4: (color online) Real-time evolution of mechanical 
quadrature X m . The thin solid (blue) line depicts an os- 
cillator with the quality factor cu m /"/ m = 10, and the thick 
solid (orange) line for u] m /^m = 100. The oscillation of the 
latter becomes stable for a longer time than the former does. 
The parameters are E/k, — 0.01, g/n = 2, u m /hi — 1, and 
Ao/io m — 1. 



simply from a time varying drive, where m (r) — 
V2g/uj m {sm(ui m T)x m - cos(a; m r)p m } and Xm(t,r) = 
— (g/uj m ) 2 sin{cj m (i — r)}. In deriving the state in (fig)) 
we have applied the decomposition in Eq. ([7]) to the sys- 
tem Hamiltonian Hs (r) in ((4]) . The oscillating functions 
<p m and Xm in the phases of Eq. (flT))) give rise to all 
harmonic components of the quantum state. Eq. (|19[) 
thus explains the cause for the dynamic quantum states 
of strongly coupled OMS. 

The damping due to the considerable rates n and 7 m 
will certainly decohere the quantum states of OMS, and 
it also determines the transient behaviors of OMS ap- 
proaching stable phases. In Figs. 2 and 3 the transient 
behavior of X c lasts for a period in the order of 1/k. 
As a comparison, we give an example of the mechanical 
quadrature X m = (x m ) evolution in Fig. 4. This quan- 
tity, which is proportional to the square of E/k, is calcu- 
lated with the averaged transformation U^(t, 0)x m U(t, 0) 
over the initial cavity vacuum state and mechanical ther- 
mal state, as well as the associate reservoir states; see 
Appendix C. One sees that it takes time in the order 
of l/7 m for mechanical quadrature X m to reach stable 
oscillation. 



C. Cavity photon number in strong coupling 
regime 

The dynamics of cavity field can be clarified further 
by the decomposition of the system-reservoir evolution 
U(t,0) = V E (t,0)V O M(t,0)V D (t,0). This decomposition 
helps to determine the contributions to a system observ- 
able by different physical factors. Below we study the 
effects of the decomposed processes on the mean cav- 
ity photon number np — (a^a). The system-reservoir 
coupling process Vb(i, 0) modifies np as it affects the 
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FIG. 5: (color online) Left: evolution of mean cavity photon number np = (a^a) in strong coupling regime; Right: mean cavity 
photon number distribution with driver detuning at a moment. The thick solid (orange) lines represent the photon number 
values including the effect of cavity noise correction in (|16(l : the dashed lines represent the photon number values under the 
approximation of the correlation function T c (t,t') = 0, where r 7^ r', in (|16[) . The cavity photon number solely generated by 
a pure drive is shown by the thin solid (blue) lines for comparison. Here the system parameters are E/k = 0.01, g/n = 2, 
i^m/«: = 2, LOraj^m = 100, and T = 0. The left frame shows the photon number evolution at the detuning point Ao = uj m , and 
the right frame depicts the photon numbers at nt = 40. 
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FIG. 6: (color online) Comparison of the cavity photon num- 
ber values np {E/k) in strong coupling regime with the corre- 
sponding predictions in the linearization approach. The dia- 
mond points represent the averaged values in the stable phase 
(see the photon number evolution pattern in Fig. 5) as calcu- 
lated in our approach, while the triangle points for the steady 
state values calculated with Eq. (6) of Ref. [Q. The pa- 
rameters of the system are chosen as g/n — 2, uim/n = 2, 
cj m /7m = 100, Ao = 0, and T = 0. 



system operators a and b in Vo(t, 0)aVs(t, 0) when act- 
ing on it; such effect exists for OMS at any temperature 
T > though it is not a major factor. The photon num- 
ber np is primarily determined by the pure drive pro- 
cess Ve(£,0) and the optomechanical coupling process 
VoM(t, 0). These main factors are described in Fig. 5 
showing the evolution of the photon number np at T = 
and its snapshot for different driver detuning at a partic- 
ular moment. Similar to the evolution of cavity quadra- 
ture, mean cavity photon number tends to a steady pe- 
riodic oscillation after a period in the order of 1/k; see 
the left frame of Fig. 5. The interference term from the 
two processes Ve(£, 0) and VoM(t, 0) contributes to the 
asymmetry of np in the blue and red detuning regime, 
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FIG. 7: (color online) Comparison between cavity photon 
numbers in bad and good cavity regime given different driver 
detuning. The left are the photon number values for ui m / k = 
0.5, and the right for uj m /K, = 3. Here the system parameters 
are E/k = 0.01, g/n = 2, uj m /j m = 100, and T = 0. These 
photon numbers are obtained at nt — 10. 



as shown in the right frame of Fig. 5. The coupling 
process VoAi(t, 0) increases np from that of a pure drive 
process Ve(£,0). For the example depicted in the fig- 
ure, it magnifies the cavity photon number by more than 
three times around the zero-detuning resonance. This is 
different from the displaced oscillator picture in [12|, LL3J , 
where the optomechanical coupling suppresses the pho- 
ton number np at the point instead. The significant de- 
viation between the cavity photon numbers calculated in 
our approach and in the linearization approach [121 Il3j 
is shown more completely in Fig. 6, which plots the pho- 
ton numbers as the function of the driving intensity E/k. 
Our approach predicts a much more significant optome- 
chanical coupling effect corresponding to the second term 
in Eq. (fT5]) . 

Contrary to the mechanical noise effect under the con- 
dition 7 m <C K, cavity noise considerably affects system 
observables in strong coupling regime. Such noise cor- 



rects the pure drive operator G(t, r) in the process of 
optomechanical coupling; see p^|) . The correction takes 
effect as long as there exists a time window in which the 
correlation function T c (t, t') of the induced cavity col- 
ored noise h c does not vanish. In Fig. 5 we compare 
the numerically calculated np with both T c (t,t / ) = 
and r c (r, t') 7^ for the drive operator in (|16l) . The 
mean cavity photon numbers in the two situations differ 
along the most of detuning range. Prominently the cavity 
photon number resonance around the red detuning point 
A = uj m and under the approximation T c (t,t') = is 
smoothed out by the cavity noise. 

Another interesting phenomenon in strong coupling 
regime is the shift of cavity photon number resonance 
from the zero detuning to the red detuned side. Given 
the same ratio g/ K, such shift will appear in the bad cav- 
ity regime with u> m /k < 1; see the example in Fig. 7. 
This could be explained with the mechanical oscillator 
displacement in proportion to g/u> m = (<?/«) X (uj m / n)~ l , 
which is larger in bad cavity so that the effective cav- 
ity resonance frequency will be shifted to a smaller value 
(note that the cavity resonance frequency is inversely pro- 
portional to cavity size). 



V. CONCLUSION 

We have studied the dynamics of OMS weakly driven in 
strong coupling regime. The damping of cavity field and 



mechanical oscillator is treated with stochastic Hamilto- 
nian for the coupling between system and reservoir [19J, 
which is also recently applied to study OMS in weak cou- 
pling regime [22|]. Different from most previous stud- 
ies, the properties of OMS described here are about the 
situations when they are genuine quantum objects, for 
example, a quantum mechanical oscillator strongly cou- 
pled to a cavity single photon, rather than their classical 
steady states plus quantum fluctuations. For the systems 
in such regimes, we illustrate the dynamical evolution of 
cavity quadrature and mean cavity photon number un- 
der weak continuous- wave drive, as well as the significant 
cavity noise effect accompanying strong optomechanical 
coupling. The method we provide is straightforwardly 
applicable to OMS driven by pulses acting for limited 
evolution time. We expect that the dynamics of general 
quantum systems with strong nonlinearity and in dissi- 
pative environment could be better understood in our 
approach as well. 
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APPENDIX 



A. Decomposition of Time-ordered Exponential 



In this appendix we provide an intuitive proof for Eq. ([7]) and a formal proof based on differential equation for Eq. 
U). Eq. ((8]) can be proved by the first method as well. 



The unitary operator Texp{— i L dT\H\(r) + ^(r))} is the infinite product of the small elements U(ti) — 
exp{—iHi(ti)5t — iH2(U)St}, where St = limjv^oo t/N and < U < t. Within the small period St, the small el- 
ement U(ti) can be decomposed into Ui(ti)Ui{bi) = Ui{U)U2{U), where Ui(tk) = exp{—iHi(tk)5t} , for any pair of 
Hi(t) and Hzit), because the cross terms due to their non-commutativity are negligible. Then we will have the 



following expression: 



Texp{-z f dr(^i(r) + #a(r))} 

Jo 



= U2{tN-l)Ul{tN-l)U2(tN-2)Ul{tN-2) • • • t^2 (*2 ) £^1 (* 2 ) ^2 (*l)^l (*l)t/ 2 (*o)^l (*o) 
= ^ 2 (*JV-l)l)l(tjV-l)^(*iV-l) 

x &a(ti\r-i)^a(*i\r-a) #i(itf-a) frj(tjv-a)#2(itf-i) 

V v ' V v ' 

V2(t,tW-2) ^(t,tj,.j) 

x U 2 {t N -i)U 2 {t N - 2 )U 2 {t N ^) Ul(t N -3) ul{t N ^)ul{t N . 2 )ul{t N ^) 
y . ' y , ' 

V2(t,t N - 3 ) v£(t,tir-3) 

X 

X U 2 (tN-l)U2(tN-2) ■ ■ ■ U 2 (tl)U2(t ) Ul{t ) U^toMih) ■ ■ ■ Ui{t N -2)Ui{t N -l) 

" v ' V v ' 

V 2 (t,0) vJ(t,o) 



cxp{-iV 2 (t,0)Hi (t a )V^ (t,0)St} 
X U2(tN-l)U2(tN-2)---U2(tl)U 2 (t ). (A-l) 

V v ' 

V 2 (t,0) 

Except for the bottom row, each row after the last equals sign in the above expression is the small element 
exp{—iV 2 (t,t k )Hi(t k )V 2 (t,t k )St}. Combining these small elements as an infinite product gives the decomposition 
formula in (|T||, 

The time-ordered exponentials Texp{— i L (It[Hi(t) + H 2 (t))} and Texp{— i L cItHi(t)} are the solutions to the 
differential equations dU/dt = —i[Hi(t) + H 2 (t))t)(t) and dVi/dt = —iHx(t)Vi(t), respectively. The initial condition 

for the differential equations is U(0) = Vi(0) = /, the identity operator. We define W(t) — V^(t)U(t). Its differential 
with respect to t reads 

_— = -V^—LvJ-U + v} =ivj H x ViV^U -iV^(H 1 +H 2 )U = -iV^H 2 ViV^U = -iV^H 2 V 1 W. (A-2) 

YAilJ LJjO UjV 

The solution of the above differential equation is W(t) — Texp{— i J Q drV^ (t)H 2 (t)Vi(t)} , implying the decomposi- 
tion formula in ©. Note that the abbreviations Vi(t,0) = Vi(i) and W(t, 0) = W(t) are used here. 

B. System Operator Transformation under System-reservoir Coupling 

We now prove the system operator transformation under Vd, the process of system-reservoir coupling in the de- 
composition of the overall system-reservoir evolution U(t,0) = Ve(i, 0)VoM(t, 0)Vd(£, 0). The general form of this 
operator transformations is V^{T 2l Ti)ciVD{T 2l Ti) , with Vd(t 2 ,ti) = Vd(^jv-i) ■ ■ ■ Vh(ti)Vh(to), an infinite product 
of the infinitesimal elements Vo(ifc) in the range t 2 > t k > t\. Because we apply the decomposition of Eq. ([7]) in the 
procedure, we will find the transformation Vd(t 2 , ti)cjVq(t2, n) as in (0), and the derivation for V^,(t 2 , ti)cjVd(t2, ti) 
is similar. Let us start with the transformation by an infinitesimal element. To do the transformation, one expands 
Voitk) out to second order in the stochastic increment [19J: 

V D {tk)ciVl(t k ) = exp{y/YidBl(tk)ci- ^/YidBi(t k )cl} di ex.p{-(y/jldBJ(tk)ci- y^dBj(*fc)c|)} 

= Ci - WTidBi{t k )c\,Ci] + 7j{VTidB}(t k )ci - y^f~dBi(t k )cl, [^fidB\(t k )c t - y/jidBi(t k )cl,Ci]] 

= {l--Jidt)c i + y/TidB i {t k ), (B-l) 

where the Ito's table 

dB\{t k )dB % {t k ) = n th dt, dBi(t k )dBj(t k ) = {n th + l)dt 

d£i(t k )d&i{t k ) = 0, dB\{t k )dB\{t k ) = (B-2) 



has been used. From this small increment we will obtain the transformation Vo(r2,Ti)ciV£)(T2,Ti) by the following 
accumulation pattern: 

Pbfo.nJcif^fo.Ti) - V D {t N -i) ■ ■ ■ Mii)VD(io)c^(*o)V^(ii) • • • ty(t N -i) 
= V D (t N -i) ■ ■ ■ V D (h)((l - pi6t)ci + y/7id&i(to))ty(ti) ■ ■ ■ V^(tiv-i) 

= Vb(tjv-i) • • • Vb(t 2 ){(l - \li5tfci + (1 - -7^i)V7^dBi(ii) + ^dBi(to)}V^(t 2 ) ■ ■ ■ P£(*jv-i) 

= (1 - \li&t) N Ci + (1 - ^^"^^(iiv-i) + • • • + (1 - ^7i**)V7i^i(*i) + VTidMto) 

= e-^- T ^l 2 c t + ^l f 2 e-^ {t '-^ )/2 Ut')dt\ (B-3) 



where dBi(t) = d L S,i{r)dr = £i(t)dt, and limjv->oo(l — ^7^)^ = e~ 7i ( T2 ~ Tl )/ 2 . Given the white noise operator £j 
satisfying [£j(t),£ [(t')] = <5(r — r'), the commutation relation [cj,cJ] = 1 is preserved under the transformation of 
Vd(t2,ti). The commutator of the term ^yi JJ 2 e~ 7i (* ~ Tl >< 2 £i(t')dt' with its conjugate does not take the form of 
delta function, so it can be regarded as a colored noise operator. 

C. System Operator Transformations under Optomechanical Coupling 

The system operator transformations under VoM{t,0) in the decomposition of the overall system-reservoir evo- 
lution U(t, 0) = Ve(i, 0)VoM{t, 0)Vd(£, 0) will be derived in this appendix. Here we use the notation VoM(t, 0) = 
Texp{— i J Q drHoM{t,T)}, with 

HoM(t,r) = -g(cos(w m T)X m (t,T) + sin(uj m T)P m (t,T))C\t,T)C(t, r) 

= -^ ro (t, ^(e-^^-^ata + e- K ( 4 - T )/ 2 aGt(t, r) + e-^-^VG^, r) + #(*, r ) G ( i; T )) . ( C -l) 

The two infinitesimal transformations for determining the transformation VQ M (t, 0)aVoM(t, 0) are 

giffoMCt.rOSr fi e -«o«(t,r,)ir = ^ + i 5 £ m ( t) ^-^-n)^) 6 + ^^(^ ri )e- K( *- Ti)/2 G(i, n)Sr (C-2) 

and 

„iH M(t,Ti)ST Q(+ T .\„ — iH M(t,Ti)&T 

= (1 + ^.g^ c (T ^ ,T,)K m (t,T ^ )G(t,T ^ )G' 1 (t,T,)<5T)G(t,T,) +i ff r c (T i ,r j )e-' c(t - Ti)/2 ^ m (i,T i )a5T, (C-3) 

where r c (ri,r,) = e~ K l n ~ Ti l/ 2 — e~ Ki - l ~ Ti } I 2 er K{ - t ~ T ^ I 2 . Also there is the commutation relation 

[K m (t, Ti), K m (t, Tj)] = [cOs(u> m Ti)X m (t, Ti) + Sm(uj m Ti)P m (t, Ti), COs(u! m Tj)X m (t, Tj) + Sm(u m Tj)P m (t, Tj)] 

— 2ie~' fm ' Ti ~ Tj '' sin u) m (Tj — Ti) = im(Ti,Tj) (C-4) 

due to the non-commutativity of X m (t, r) and P m (t, t). 

One of the advantages in our approach of combined system-reservoir unitary evolution is the availability of the 
property U^(t,0)f(c u c\)U(t,0) = f(U^(t,0)diU(t,0),U^(t,0)c\U(t,0)) of a unitary transformation. This property 
enables one to obtain a formal expression for the transformation Vq M (£, 0)aVoM(t,0) as follows: 

iH M(t,TN-l)ST _ _ _ iH M(t,Tl)ST iHoM(t,T )ST •> e — iHoM{t,T )ST -iH M(t,Tl)ST _ _ _ -iH M(t,TN-l)ST 

= (1 + igJkg(t, to^-^-^St) (1 + igJkg(t, t^'^'^St) ■ ■ ■ (l + tfl^(t, Tjv-Oe-"**-™- 1 ^ 
+ ig(l +igkg{t,To)e- K{ - t - T ^8T) ••• (l + ^.g^(i,T W _ 2 ) e - K ( t - T "- 2 ) ( 5T)^(i,T W _ 1 )Gc(t,T J v_ 1 ) ( 5T 
+ ••• + *. 9 (l + i3 i^(t,To)e- K ( t -^T) e - K ^^ 

= (Texp{-z 5 fdre-^-^K^r)})^ 
Jo 

+ i s / dr (Tcxp{-* 3 f du e-^-^Kgit^e-^-^^Kgit^Gc^T), (C-5) 
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where K%(t,Ti) = Vl M {t,Ti)k m {t,Ti)VoM{t,n)-, and G c (t,n) = v£ M (t,Ti)G(t,n)VoM{t,Ti). 

To apply (|C-5|) to numerical calculations, one should find the proper forms of K^(t, r) and Gc(t, r). The modified 
operator K^it, n) is obtained by the successive infinitesimal unitary operations (from tj+i to t) on K m (t,Ti), which 
comes from the infinitesimal transformation e lH oM(t,Ti)Sr » e -%H M(t,n)ST a ^ ^ e momen t Ti - see the following: 

J{ c (f T .) = p?HoM(t,TN-\)&T ___ iH M{t,T i+ i)Sr j£ i± \ —iHoM(t,n+i)ST _ _ .—iHoM(t,TN-if>r 

= k m (t,Ti) + g f dTm(r,Ti)V^ M (t,T)CHt,T)C(t,T)VoM(t,T). (C-6) 

The above equation is expanded to 

kg{t,n) = K m (t,Ti)+g f dTe- K ^-^m(T,Ti)V^ M (t,T)a^aVoM{t,T) 

+ 9 I At e - K ( 4 - T )/ 2 TO(^^)^Af(i^)(aGt(t,r)+atG(i,T))V> OM (t,r) 

+ g f dr m(T,n)ty M (t,T)Gi(t,T)G(t,T)V M(t,T). (C-7) 



The term carrying the decayed correlation function e K (* T 'm(T, Ti) in the integrand can be neglected. 
Next, based on the infinitesimal transformation (|C-3[I . one will also find 



G c {t,n) = G(t,n)+ig dr r c (r,rO^ M (i,r)^ m (i, r)f OA f (i,T)V^ M (i,T)G(i,T)T/ OM (i,T) 

J Ti 

+ ig f dre-"(^)/ 2 r c (r,r i )^(t,r)it n (t,r)t^M(*,r)^(t,r)6 1 PbM(*,r). (C-8) 



The last term with the decayed correlation function e K '* T '/ 2 r c (-r, Tj) in the integrand can be well neglected, and 
the remaining terms will be expanded by the iteration of the above to 



G c (t,n) = G(t,n)+ig f dTT c (T,Ti)Kg(t,T)G(t,T) 

ft pt 



+ {igf / drr c (T,n)Kg(t,r) / drT c (r' »*£(*, /)<?(«, /) 

+ ( ?5 ) 3 [* dTT c ( T ,n)kC(t,T) f dT'r c (r',T)kg(t,T') [ t dT"r c (r",T')k^(t,T")G(t,T") + 



(C-9) 

In principle, the negligence of the integrals carrying the decayed correlation functions for Eqs. (j(J-7[) and (|CJ-8[) is 
the only approximation made in our procedure. For a point n outside the vicinity of the end time point t, the decay 
factor e _K (*~ T )/ 2 completely damps the correlation function T c (t, Ti) inside the correlation time window around it to 
zero, so there is no contribution from the neglected integral. The correlation function m(r, Tj) defined in (IC-4|) takes 
the form of oscillating function for large quality factor u) m /'y m 3> 1, suppressing the integral of the positive term 
Vq M (P, T)a)aVoM{t,T) even without the damping factor e _K '* _T ^. 

More approximations could be made to simplify the above expressions. First, the summation of multi-folded 
integrals in (|C-9p can be conveniently estimated given a short correlation time window of the colored cavity noise, in 
which the drive operator G(t, r) changes slowly. In this way such summation can be approximated by a time-ordered 
exponential. Second, for a weak drive of E/k <C 1, the operator k^(t,Ti) can be simply approximated by k m (t,Ti). 
Putting these approximations together, one will have the following closed form of the transformation 

v£ M (t,0)aV O M(t,Q) 

_ e tg p dre-^-T) S i a (u m r)P m (t,r) e tg p dre-^-T) cos (^ m T)X m {t,r) e -ig 2 Q{t) ^ 

rt 

+ iq c / jTe - i 9 2 ®( T ) e i 9 So due ~" (t ~ u} sin(u] m u)P m (t,u) e ig fo due - "** - "' cos(u m u)X m (t,u) 

X e~ K( ' t ~ T ^ 2 k m (t T \S e i sf*'i^ T o(u,T)am(u m u)P m (t,u) e igf*dur c (u,T)cos(u m u)X nt (t,u) e -ig 2 e'(T)Qfj. nj (C-10) 



11 



for a single-photon weak driver. Here we also have the extra phases 



e(r)=2/ due- K ^- u hm(uj m u) dve -^t-v)- lm {u-v)/2 CO s( Wm «), 
Jo Jo 



e'(r) = 2/ duT c (u,r)sia(uj m u) dv T c (v,T)e^ miu - v)/2 cos(uj m v) 



after factorizing the time-ordered exponentials involving X m (t, u) and P m (t, u). 
Moreover, the pure driving operation on cavity operator is simply found as 

V£{t,0)aV E (t,0) = Texp{£ f dT(ae _iA ° T e-t (t - r) - h.c.)} aTexpjfi f dT(aV AoT e-t ( *- T) - /i.e.)} 

Jo Jo 

= a + E [ dr e * A ° r e-*( t - T > = a + £>i(t). (C-ll) 

Jo 

Finally, the transformations of the mechanical oscillator operators W(t, 0)x m U(t,0) and W(t,0)p m U(t,0), where 
U(t,0) — Ve(£, 0)Vc»M(t, 0)Vb(t, 0), can be found in a similar way to (|C-6[) . Under the same approximation as for 
(|C-10[) , their averages over the initial state of cavity vacuum, thermal state of mechanical oscillator, as well as the 
cavity vacuum reservoir and mechanical thermal reservoir state, read 



(U\t,0)x m U(t,0)) = -V2g / dr e-^ t -^' 2 ^{u Jm T)(e~^ t ^l 2 D* 1 {T) + D* 2 (r))(e-^ t -^/ 2 D 1 {r) + D 2 {t)), 

Jo 

(U*(t,0)p m U(t,0)) = V2g f dre-^ t -^l 2 cos{u: m T){e-^ t -^l 2 Dl{T)+D* 2 {T)){e- K{ - t -^l 2 D l {T)+D 2 {T)). (C-12) 



These averages are used to study mechanical oscillator dynamics under general optomechanical coupling. 
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